#include <limits>

#include <QtCore/QFile>
#include <QtCore/QTextStream>
#include <QtCore/QStringList>

#include "glwidget.h"
#include "geometry/2d/point.h"
#include "fem/trianglenet.h"
#include "fem/frontaltriangulation.h"
#include "fem/delaunaytriangulation.h"
#include "approximation/interpolation.h"

//================================================================================================================================
//Отрисовка OpenGL. Вызывается автоматически при отрисовке виджета.
//================================================================================================================================
void GLWidget::paintGL()
{
    glClear(GL_COLOR_BUFFER_BIT);

    //****************************************************************************************************************************
    //Устанавливаем масштаб.
    //****************************************************************************************************************************
    glMatrixMode(GL_PROJECTION);
    glLoadIdentity();
    gluOrtho2D(m_Scene.left, m_Scene.right, m_Scene.bottom, m_Scene.top);

    //****************************************************************************************************************************
    //Рисуем сетку.
    //****************************************************************************************************************************
    glCallList(m_NetListID);

    //****************************************************************************************************************************
    //Рисуем контур тигля.
    //****************************************************************************************************************************
    glCallList(m_AreaListID);

    glFinish();
}

//================================================================================================================================
//================================================================================================================================
void GLWidget::paintMesh(int method)
{
    if(method == 0)
    {
        fem::FrontalTriangulation triangulation(QStringList() << "T" << "Vf" << "Psi" << "Omega");

        int base_points[15];
        base_points[0] =  triangulation.addBasisPoint(geometry::_2d::Point(0.000, 0.000));
        base_points[1] =  triangulation.addBasisPoint(geometry::_2d::Point(0.055, 0.000));
        base_points[2] =  triangulation.addBasisPoint(geometry::_2d::Point(0.000, 0.005));
        base_points[3] =  triangulation.addBasisPoint(geometry::_2d::Point(0.050, 0.005));
        base_points[4] =  triangulation.addBasisPoint(geometry::_2d::Point(0.000, 0.055));
        base_points[5] =  triangulation.addBasisPoint(geometry::_2d::Point(0.050, 0.055));
        base_points[6] =  triangulation.addBasisPoint(geometry::_2d::Point(0.055, 0.055));
        base_points[7] =  triangulation.addBasisPoint(geometry::_2d::Point(0.000, 0.060));
        base_points[8] =  triangulation.addBasisPoint(geometry::_2d::Point(0.050, 0.060));
        base_points[9] =  triangulation.addBasisPoint(geometry::_2d::Point(0.010, 0.065));
        base_points[10] = triangulation.addBasisPoint(geometry::_2d::Point(0.055, 0.065));
        base_points[11] = triangulation.addBasisPoint(geometry::_2d::Point(0.000, 0.115));
        base_points[12] = triangulation.addBasisPoint(geometry::_2d::Point(0.010, 0.115));
        base_points[13] = triangulation.addBasisPoint(geometry::_2d::Point(0.000, 0.006));
        base_points[14] = triangulation.addBasisPoint(geometry::_2d::Point(0.015, 0.005));

        int base_segments[19];
        base_segments[0] =  triangulation.addBorderSegment(base_points[0], base_points[1], 100);
        base_segments[1] =  triangulation.addBorderSegment(base_points[0], base_points[2], 6);
        base_segments[2] =  triangulation.addBorderSegment(base_points[1], base_points[3], 6);
        base_segments[3] =  triangulation.addBorderSegment(base_points[2], base_points[3], 80, 0.98);
        base_segments[4] =  triangulation.addBorderSegment(base_points[2], base_points[4], 80, 1.02);
        base_segments[6] =  triangulation.addBorderSegment(base_points[1], base_points[6], 100);
        base_segments[8] =  triangulation.addBorderSegment(base_points[5], base_points[6], 6);
        base_segments[9] =  triangulation.addBorderSegment(base_points[4], base_points[7], 6);
        base_segments[10] = triangulation.addBorderSegment(base_points[5], base_points[8], 6);
        base_segments[11] = triangulation.addBorderSegment(base_points[6], base_points[10],10);
        base_segments[12] = triangulation.addBorderSegment(base_points[7], base_points[8], 60);
        base_segments[13] = triangulation.addBorderSegment(base_points[8], base_points[10],6);
        base_segments[14] = triangulation.addBorderSegment(base_points[7], base_points[11],60);
        base_segments[15] = triangulation.addBorderSegment(base_points[7], base_points[9], 8);
        base_segments[16] = triangulation.addBorderSegment(base_points[9], base_points[10],60);
        base_segments[17] = triangulation.addBorderSegment(base_points[9], base_points[12],60);
        base_segments[18] = triangulation.addBorderSegment(base_points[11],base_points[12],8);

        triangulation.addQuadrangleSubdomain(base_segments[13], base_segments[11], base_segments[8],  base_segments[10],0);
        triangulation.addQuadrangleSubdomain(base_segments[16], base_segments[13], base_segments[12], base_segments[15],0);
        triangulation.addQuadrangleSubdomain(base_segments[18], base_segments[17], base_segments[15], base_segments[14],0);

        int add_points[2];
        add_points[0] = triangulation.addBasisPoint(geometry::_2d::Point(0.000, 0.006));
        add_points[1] = triangulation.addBasisPoint(geometry::_2d::Point(0.025, 0.005));

        double r_first(0.000), r_last(0.025);
        double z_first(0.006), z_last(0.005);

        QList<double> rs, zs;
        triangulation.getSegmentPoints(3, rs, zs);

        int j = 0;
        while(rs[j] < 0.025)
            ++j;
        while(j < rs.size())
            rs.removeLast();

        QList<int> front_plus_indexes;
        QList<double> front_plus_interpolated = approximation::Interpolation::interpolate(QList<double>() << 0.000 << 0.025, QList<double>() << 0.006 << 0.005, rs, approximation::Interpolation::Linear);
        front_plus_indexes.append(add_points[0]);
        for(int i = 1; i < front_plus_interpolated.size() - 1; i++)
            front_plus_indexes.append(triangulation.addBasisPoint(geometry::_2d::Point(rs[i],front_plus_interpolated[i])));
        front_plus_indexes.append(add_points[1]);

        int add_segments[7];
        add_segments[0] = triangulation.addBorderSegment(base_points[2],add_points[1],front_plus_indexes.size() - 2, 0.98);
        add_segments[1] = triangulation.addBorderSegment(add_points[1],base_points[3],int(80*(0.050 - r_last)/0.05), 0.98);
        add_segments[2] = triangulation.addBorderSegment(base_points[2],add_points[0],int(80*(z_first - 0.005)/0.05), 1.02);
        add_segments[3] = triangulation.addBorderSegment(add_points[0],base_points[4],int(80*(0.055 - z_first)/0.05), 1.02);
        add_segments[4] = triangulation.addBorderSegment(fem::BorderSegment(front_plus_indexes, 0.98));
        add_segments[5] = triangulation.concatinateBorderSegments(add_segments[0],add_segments[1]);
        add_segments[6] = triangulation.concatinateBorderSegments(add_segments[4],add_segments[1]);

        base_segments[5] = triangulation.addBorderSegment(base_points[3],base_points[5],int(80*(z_first - 0.005)/0.05) + int(80*(0.055 - z_first)/0.05), 1.02);
        base_segments[7] = triangulation.addBorderSegment(base_points[4],base_points[5],front_plus_indexes.size() - 2 + int(80*(0.050 - r_last)/0.05) + 1, 0.98);

        triangulation.addQuadrangleSubdomain(base_segments[8], base_segments[6], base_segments[2], base_segments[5], 0);
        triangulation.addQuadrangleSubdomain(base_segments[12], base_segments[10], base_segments[7], base_segments[9], 3);
        triangulation.addQuadrangleSubdomain(add_segments[5], base_segments[2], base_segments[0], base_segments[1], 0);
        triangulation.addQuadrangleSubdomain(base_segments[7], base_segments[5], add_segments[6], add_segments[3], 1);
        triangulation.addTriangleSubdomain(add_segments[0], add_segments[4], add_segments[2], 4);

        fem::TriangleNet net = triangulation.triangulate();
        double x_min(DBL_MAX), x_max(DBL_MIN), y_min(DBL_MAX), y_max(DBL_MIN);
        for(int i = 0; i < net.numOfPoints(); ++i)
        {
            if(net.point(i).x() < x_min) x_min = net.point(i).x();
            if(net.point(i).x() > x_max) x_max = net.point(i).x();
            if(net.point(i).y() < y_min) y_min = net.point(i).y();
            if(net.point(i).y() > y_max) y_max = net.point(i).y();
        }

        m_Scene.left = x_min*0.98;
        m_Scene.bottom = y_min*0.98;
        m_Scene.right = x_min + qMax(x_max - x_min, y_max - y_min)*1.03;
        m_Scene.top = y_min + qMax(x_max - x_min, y_max - y_min)*1.03;

        m_Scene.dtop = m_Scene.top;
        m_Scene.dleft = m_Scene.left;
        m_Scene.dright = m_Scene.right;
        m_Scene.dbottom = m_Scene.bottom;

        //****************************************************************************************************************************
        //****************************************************************************************************************************
        glDeleteLists(m_NetListID, 1);
        glDeleteLists(m_AreaListID, 1);
        m_NetListID = glGenLists(1);
        m_AreaListID = glGenLists(1);

        glNewList(m_NetListID, GL_COMPILE);
            qglColor(Qt::black);
            for(int i = 0; i < net.numOfElements(); ++i)
            {
                glBegin(GL_LINE_LOOP);
                for(int j = 0; j < 3; ++j)
                    glVertex2f(net.point(i,j).x(), net.point(i,j).y());
                glEnd();
            }
        glEndList();

        glNewList(m_AreaListID, GL_COMPILE);
            qglColor(Qt::black);
            glLineWidth(2);
            glBegin(GL_LINES);
                glVertex2d(0.000, 0.005);
                glVertex2d(0.050, 0.005);
                glVertex2d(0.050, 0.005);
                glVertex2d(0.050, 0.060);
                glVertex2d(0.000, 0.055);
                glVertex2d(0.050, 0.055);
                glVertex2d(0.000, 0.060);
                glVertex2d(0.050, 0.060);
                glVertex2d(0.000, 0.000);
                glVertex2d(0.055, 0.000);
                glVertex2d(0.055, 0.000);
                glVertex2d(0.055, 0.065);
                glVertex2d(0.055, 0.065);
                glVertex2d(0.010, 0.065);
                glVertex2d(0.010, 0.065);
                glVertex2d(0.010, 0.115);
                glVertex2d(0.010, 0.115);
                glVertex2d(0.000, 0.115);
                glVertex2d(0.000, 0.115);
                glVertex2d(0.000, 0.000);
            glEnd();
            glLineWidth(1);
        glEndList();
    }
    //****************************************************************************************************************************
    //****************************************************************************************************************************
    if(method == 1)
    {
        fem::DelaunayTriangulation triangulator(QStringList("T"));

        int points[4];
        points[0] = triangulator.addPoint(geometry::_2d::Point(0.000, 0.000));
        points[1] = triangulator.addPoint(geometry::_2d::Point(1.000, 0.000));
        points[2] = triangulator.addPoint(geometry::_2d::Point(0.000, 1.000));
        points[3] = triangulator.addPoint(geometry::_2d::Point(1.000, 1.000));

        int segments[4];
        segments[0] = triangulator.addBorderSegment(points[0], points[1], 25);
        segments[1] = triangulator.addBorderSegment(points[0], points[2], 25);
        segments[2] = triangulator.addBorderSegment(points[1], points[3], 25);
        segments[3] = triangulator.addBorderSegment(points[2], points[3], 25);

        int superstructures[1];
        superstructures[0] = triangulator.addSuperstructure(QList<int>() << points[0] << points[1] << points[2] << points[3], 0);

        for(int i = 0; i < 300000; i++)
            triangulator.addPoint(superstructures[0], geometry::_2d::Point(double(qrand())/RAND_MAX, double(qrand())/RAND_MAX));

        fem::TriangleNet net = triangulator.triangulate();

        double x_min(DBL_MAX), x_max(DBL_MIN), y_min(DBL_MAX), y_max(DBL_MIN);
        for(int i = 0; i < net.numOfPoints(); ++i)
        {
            if(net.point(i).x() < x_min) x_min = net.point(i).x();
            if(net.point(i).x() > x_max) x_max = net.point(i).x();
            if(net.point(i).y() < y_min) y_min = net.point(i).y();
            if(net.point(i).y() > y_max) y_max = net.point(i).y();
        }

        m_Scene.left = x_min*0.98;
        m_Scene.bottom = y_min*0.98;
        m_Scene.right = x_min + qMax(x_max - x_min, y_max - y_min)*1.03;
        m_Scene.top = y_min + qMax(x_max - x_min, y_max - y_min)*1.03;

        m_Scene.dtop = m_Scene.top;
        m_Scene.dleft = m_Scene.left;
        m_Scene.dright = m_Scene.right;
        m_Scene.dbottom = m_Scene.bottom;

        //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
        //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
        glDeleteLists(m_NetListID, 1);
        glDeleteLists(m_AreaListID, 1);
        m_NetListID = glGenLists(1);
        m_AreaListID = glGenLists(1);

        glNewList(m_NetListID, GL_COMPILE);
            qglColor(Qt::black);
            for(int i = 0; i < net.numOfElements(); ++i)
            {
                glBegin(GL_LINE_LOOP);
                for(int j = 0; j < 3; ++j)
                    glVertex2f(net.point(i,j).x(), net.point(i,j).y());
                glEnd();
            }
        glEndList();

        glNewList(m_AreaListID, GL_COMPILE);
            qglColor(Qt::black);
            glLineWidth(2);
            glBegin(GL_LINES);
                glVertex2d(0.000, 0.000);
                glVertex2d(1.000, 0.000);
                glVertex2d(1.000, 1.000);
                glVertex2d(0.000, 1.000);
                glVertex2d(0.000, 0.000);
            glEnd();
            glLineWidth(1);
        glEndList();
    }

    //****************************************************************************************************************************
    //****************************************************************************************************************************
    update();
}

//================================================================================================================================
//================================================================================================================================
void GLWidget::paintMesh(QString filename)
{
    int n(0);
    QFile file(filename);
    QTextStream istream(&file);
    fem::FrontalTriangulation triangulator(QStringList("T"));

    QList<int> base_points, add_points;
    QList<int> base_segments, add_segments;

    file.open(QFile::ReadOnly | QFile::Text);
    istream >> n;
    for(int i = 0; i < n; ++i)
    {
        double x, y;
        istream >> x >> y;
        base_points.append(triangulator.addBasisPoint(geometry::_2d::Point(x, y)));
    }
    istream >> n;
    for(int i = 0; i < n; ++i)
    {
        double factor;
        int begin, end, num_of_points;
        istream >> begin >> end >> num_of_points >> factor;
        base_segments.append(triangulator.addBorderSegment(begin, end, num_of_points, factor));
    }
    istream >> n;
    for(int i = 0; i < n; ++i)
    {
        int segment_1, segment_2;
        istream >> segment_1 >> segment_2;
        triangulator.concatinateBorderSegments(segment_1, segment_2);
    }
    istream >> n;
    for(int i = 0; i < n; ++i)
    {
        int segments[3], material;
        for(int j = 0; j < 3; ++j)
            istream >> segments[j];
        istream >> material;
        triangulator.addTriangleSubdomain(segments[0], segments[1], segments[2], material);
    }
    istream >> n;
    for(int i = 0; i < n; ++i)
    {
        int segments[4], material;
        for(int j = 0; j < 4; ++j)
            istream >> segments[j];
        istream >> material;
        triangulator.addQuadrangleSubdomain(segments[0], segments[1], segments[2], segments[3], material);
    }
    file.close();

    add_points.append(triangulator.addBasisPoint(geometry::_2d::Point(0.000, 0.006)));
    add_points.append(triangulator.addBasisPoint(geometry::_2d::Point(0.025, 0.005)));

    QList<double> rs, zs;
    triangulator.getSegmentPoints(3, rs, zs);

    int j = 0;
    while(rs[j] < 0.025)
        ++j;
    while(j < rs.size())
        rs.removeLast();

    QList<int> front_plus_indexes;
    QList<double> front_plus_interpolated = approximation::Interpolation::interpolate(QList<double>() << 0.000 << 0.025, QList<double>() << 0.006 << 0.005, rs, approximation::Interpolation::Linear);
    front_plus_indexes.append(add_points[0]);
    for(int i = 1; i < front_plus_interpolated.size() - 1; i++)
        front_plus_indexes.append(triangulator.addBasisPoint(geometry::_2d::Point(rs[i],front_plus_interpolated[i])));
    front_plus_indexes.append(add_points[1]);

    add_segments.append(triangulator.addBorderSegment(base_points[2],add_points[1],front_plus_indexes.size() - 2, 0.98));
    add_segments.append(triangulator.addBorderSegment(add_points[1],base_points[3],(int)(80*(0.050 - 0.025)/0.05), 0.98));
    add_segments.append(triangulator.addBorderSegment(base_points[2],add_points[0],(int)(80*(0.006 - 0.005)/0.05)));
    add_segments.append(triangulator.addBorderSegment(add_points[0],base_points[4],(int)(80*(0.055 - 0.006)/0.05)));
    add_segments.append(triangulator.addBorderSegment(fem::BorderSegment(front_plus_indexes)));
    add_segments.append(triangulator.concatinateBorderSegments(add_segments[0],add_segments[1]));
    add_segments.append(triangulator.concatinateBorderSegments(add_segments[4],add_segments[1]));

    triangulator.addQuadrangleSubdomain(add_segments[5], base_segments[2], base_segments[0], base_segments[1], 0);
    triangulator.addQuadrangleSubdomain(base_segments[7], base_segments[5], add_segments[6], add_segments[3], 1);
    triangulator.addTriangleSubdomain(add_segments[0], add_segments[4], add_segments[2], 4);

    fem::TriangleNet net = triangulator.triangulate();
    double x_min(DBL_MAX), x_max(DBL_MIN), y_min(DBL_MAX), y_max(DBL_MIN);
    for(int i = 0; i < net.numOfPoints(); ++i)
    {
        if(net.point(i).x() < x_min) x_min = net.point(i).x();
        if(net.point(i).x() > x_max) x_max = net.point(i).x();
        if(net.point(i).y() < y_min) y_min = net.point(i).y();
        if(net.point(i).y() > y_max) y_max = net.point(i).y();
    }

    m_Scene.left = x_min*0.98;
    m_Scene.bottom = y_min*0.98;
    m_Scene.right = x_min + qMax(x_max - x_min, y_max - y_min)*1.03;
    m_Scene.top = y_min + qMax(x_max - x_min, y_max - y_min)*1.03;

    m_Scene.dtop = m_Scene.top;
    m_Scene.dleft = m_Scene.left;
    m_Scene.dright = m_Scene.right;
    m_Scene.dbottom = m_Scene.bottom;

    //****************************************************************************************************************************
    //****************************************************************************************************************************
    glDeleteLists(m_NetListID, 1);
    glDeleteLists(m_AreaListID, 1);
    m_NetListID = glGenLists(1);
    m_AreaListID = glGenLists(1);

    glNewList(m_NetListID, GL_COMPILE);
        qglColor(Qt::black);
        for(int i = 0; i < net.numOfElements(); ++i)
        {
            glBegin(GL_LINE_LOOP);
            for(int j = 0; j < 3; ++j)
                glVertex2f(net.point(i,j).x(), net.point(i,j).y());
            glEnd();
        }
    glEndList();

    glNewList(m_AreaListID, GL_COMPILE);
        qglColor(Qt::black);
        glLineWidth(2);
        glBegin(GL_LINES);
            glVertex2d(0.000, 0.005);
            glVertex2d(0.050, 0.005);
            glVertex2d(0.050, 0.005);
            glVertex2d(0.050, 0.060);
            glVertex2d(0.000, 0.055);
            glVertex2d(0.050, 0.055);
            glVertex2d(0.000, 0.060);
            glVertex2d(0.050, 0.060);
            glVertex2d(0.000, 0.000);
            glVertex2d(0.055, 0.000);
            glVertex2d(0.055, 0.000);
            glVertex2d(0.055, 0.065);
            glVertex2d(0.055, 0.065);
            glVertex2d(0.010, 0.065);
            glVertex2d(0.010, 0.065);
            glVertex2d(0.010, 0.115);
            glVertex2d(0.010, 0.115);
            glVertex2d(0.000, 0.115);
            glVertex2d(0.000, 0.115);
            glVertex2d(0.000, 0.000);
        glEnd();
        glLineWidth(1);
    glEndList();

    //****************************************************************************************************************************
    //****************************************************************************************************************************
    update();
}
